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SUMMARY 


Finite-difference solutions were developed for treating one-dimensional 
transient diffusion in single-, two-, and three-phase binary alloy systems. 

These solutions are applicable for planar, cylindrical, or spherical geometries 
with any diffusion-zone size and any continuous (within each phase) variation 
of the diffusion coefficient with concentration. Special techniques were 
included to account for differences in molal volumes, initiation and growth of 
an intermediate phase (three-phase system), disappearance of a phase (two- and 
three-phase systems), and the presence of an initial composition profile in the 
specimen. In each analysis, an effort was made to achieve good accuracy while 
minimizing computation time. A major improvement in solution accuracy was 
achieved in the two-phase analysis by employing a mass-conservation criterion to 
establish the location of the interface rather than the conventional interface- 
flux-balance criterion. In the three-phase analysis, computation time was mini- 
mized without sacrificing solution accuracy by treating the three-phase problem 
as a two-phase problem when the thickness of the intermediate phase was less 
than a preset small value. The computation time was also minimized, without a 
significant loss of accuracy, by increasing the grid size for longer diffusion 
times when the concentration in each phase was relatively uniform. For the 
single-phase analysis, this was accomplished by reducing the number of grid 
stations with diffusion time. For the two- and three-phase analyses, the grid 
stations were initially concentrated near the interface, and the zone of analy- 
sis was expanded with diffusion time. Comparisons with other solutions were 
made where possible. 


INTRODUCTION 

Solutions of the diffusion equation have been reported (refs. 1 to 7) for 
several different initial and boundary conditions. Most of these, however, are 
restricted to one geometry, applicable only for infinite or semi-infinite sys- 
tems, or require that the diffusion coefficient (D) be constant. Heckel and 
coworkers (refs. 1 to 3) developed finite-difference solutions for two- and 
three-phase binary alloy systems capable of treating different geometries and 
any diffusion-zone size, but they assumed D to be a constant. Unnam and 
Houska (refs. ^ and 5) have reported iterative solutions for single- and two- 
phase systems allowing for variation of D with concentration, but these solu- 
tions were restricted to planar geometry. Finite-difference solutions for 
single-phase (ref. 6) and two-phase (ref. 7) binary alloy systems allowing for 
a concentration-dependent D have been reported. However, these solutions 
were developed only for cylindrical geometry. The computer codes for even 
these restricted solutions have not, in general, been made available. There- 
fore, a need exists for more general solutions and their computer codes. 

The purpose of this report is to describe general finite-difference solu- 
tions for single-, two-, and three-phase binary alloy systems which can treat 


one-dimensional diffusion in planar, cylindrical, or spherical geometry speci- 
mens with any diffusion-zone size and any concentration dependence of the dif- 
fusion coefficient. General computer codes were developed to perform these 
analyses and are available from COSMIC (refs. 8, 9, and 10). Also available 
with each code are documentation describing the code, a program listing, and 
a sample case illustrating typical input and output information. The special 
features of each analysis are discussed in this report, and example cases are 
presented to illustrate the validity of the solutions. 


SYMBOLS 

elements which constitute binary alloy system 
composition, atomic fraction of B 

average composition for sample, atomic fraction of B 

initial concentration of 8 -phase, assumed to be unity in this study, 
atomic fraction of B 

Cq initial concentration of a-phase, assumed to be zero in this study, 

atomic fraction of B 

D chemical diffusion coefficient, m^/s 

lA chemical diffusion coefficient of pure A, m^/s 

D® chemical diffusion coefficient of pure B, m^/s 

D normalized diffusion coefficient, D/Dmax 

D mav maximum diffusion coefficient at temperature of interest for entire 

composition range in sample, m^/s 

1 number of grid stations in 8 -phase in two-phase system including grid 

station at 8ot interface 

Kjj normalized distance in the yphassi s®® equation (26) 

L total thickness or diameter of specimen, m 

% initial thickness or diameter of B-rich region, m 

m parameter whose value depends on type of geometry: 0 for planar, 

1 for cylindrical, and 2 for spherical 

N total number of grid stations 

n parameter denoting finite-difference grid station 

R normalized distance, r/(L/2) 
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A,B 

C 

C 

C 


r distance from center of B-rich region; a subscript n denotes the 

distance to nth finite-difference grid station, m 

T normalized time, P mq yb/ (L/2)^ 

t time , s 

a, 3 designations for terminal phases rich in A and B, respectively 

Y designation for intermediate phase 

difference between solubility limits of 3- and a-phases, C^qj - 
(other double subscripts have similar meaning), atomic fraction 
of B 

©D temperature of interest, K 

5/2 distance from center of 3-phase to 3a interface in two-phase system, 

positive superscript designates a side of interface and negative 
superscript designates 3 side of interface, m 

5 1/2 distance from center of 3 -phase to 3 y interface in three-phase sys- 

tem, positive superscript designates y side of interface and nega- 
tive superscript designates 3 side of interface, m 

C 2/2 distance from center of 3 -phase to y“ interface in three-phase sys- 

tem, positive superscript designates a side of interface and nega- 
tive superscript designates y side of interface, m 

Superscripts: 

j index of time increment 

a A-rich phase 

3 B-rich phase 

Y intermediate phase 
Subscripts: 

I grid-station number in 3-phase at 3a interface of two-phase system 

I-] grid-station number in 3-phase at 3 y interface of three-phase system 

I 2 grid-station number in Y~Phase at ya. interface in three-phase system 

N Nth grid station of finite-difference grid 

n index of grid station of finite-difference grid 
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Double subscripts o3 , 3a, ay, etc., designate concentration or diffusion 

coefficient at a particular interface. The first letter indicates the phase 
being referenced and the two letters together indicate the interface being con- 
sidered. For example, Qjg denotes the concentration in the a-phase at the 
a3 interface. 


THEORETICAL DEVELOPMENT 
Single-Phase Systems 

Diffusion between pure metals A and B which belong to a single-phase binary 
alloy system will result in a continuous composition variation of the form shown 
in figure 1(a). Diffusion in this system can be described by one-dimensional 
Pick's second law 



where C is the atomic fraction of B, r is the distance from the center of 
the B-rich region, D is the concentration-dependent diffusion coefficient, and 
m = 0, 1, or 2 for planar, cylindrical, or spherical geometries, respectively. 
The initial and boundary conditions to be satisfied are 

C = C (0 ^ r < il/2; t = 0) 

C = Co (S</2 < r < L/2; t = 0) 

^ =0 (r = 0 and L/2; t ^ 0) 

9r 

where C and Cq are the initial compositions of the B-rich and A-rich alloys 
joined to form the diffusion couple, A is the initial thickness or diameter of 
the B-rich region, and L is the total thickness or diameter of the specimen. 

The finite-difference expression (explicit form, second-order central dif- 
ference) for Pick's second law is given by 


pj 

^n ” "^n 
At 


= dJ 


m 


^ cj_i ^ 

\ 2 Ar / \ (Ar)2 / 


fr^ 1 - 1 

*^n+1 ^n-1 

( 2 Ar / 


'nJ nJ ' 
%+1 ~ Pn -1 

{ 2 Ar 


For r=0 (n=1), the finite-difference expression is given by 

j+1 j / j J'' 

= 2(m + 1 )DjF 2 .Jl_C _1 
At 1 \ (Ar)2 


( 2 ) 


(3) 
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To obtain this expression, a second-order forward difference was used, and the 
singularity in the differential equation (eq. (1) for m = 1 or 2) at r = 0 
was taken into account (ref. 11). At r = L/2 (n = N), equation (2) reduces 
to 

oJ -i (J 

- Cn = ^N-1 - Cn (4) 

it (ir)2 / 

The j superscripts designate time increments and the n subscripts designate 
grid stations in the diffusion couple. 

Finite-difference calculations were made by choosing the grid spacing small 
enough to give convergent solutions of acceptable accuracy. Because explicit 
formulation was used in the development of the equations, the time increment 
At was selected according to the stability requirement. 

At i 0.25 (5) 

^max 

where Djuax is the maximum diffusion coefficient for the system at the tempera- 
ture of interest. 

A computer code was developed to solve equations (2) to (5) for any binary 
alloy system which exhibits complete solid solubility. The FORTRAN IV code, 
documentation describing the code, and a sample case illustrating typical input 
and output information are available from COSMIC (ref. 8). The special features 
of the analysis included in this code are discussed in a subsequent section. 


Two-Phase Systems 


Diffusion between pure metals which form a two-phase binary alloy system 
will result in a composition variation between the two metals of the form shown 
in figure 1(b). The concentrations at the a8 interface correspond to the 
equilibrium solubility limits of the 3-phase (C^qj) and the a-phase (C^jjg) at 
the temperature 9 d of interest. 


Diffusion in a two-phase system can be described by two partial differen- 
tial equations (Fick's second law for each phase) and an interface-flux-balance 
equation which describes motion of the interface. The interface-flux-balance 
equation is given by 




- Ife ('oc*' 

V=C“/2 


(6) 


where ^/2 is the location of the a3 interface, with superscripts + and - 
designating the a and 3 sides of the interface, respectively. The initial 
and boundary conditions to be satisfied are 


C = C* 
C = 


(0 1 r < il/2; t = 0) 
(il/2 < r % L/2; t = 0) 
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^ = 0 
ar 


(r = 0 and L/2; t ^ 0) 


C = Cfl 


(r = ^'/2) 


C = C, 


(r = ^■"/2) 


In general, the thicknesses of the a- and 3-phases change with diffusion 
time. To account for these dimensional changes, the variable-grid technique 
developed by Murray and Landis (ref. 12) was used. The number of stations in 
each phase was held constant, and the grid sizes Ar® and Ar^ changed with 
time. The 3-phase was divided into I - 1 equally sized space increments of 
thickness ArP . The a-phase was divided into N - I - 1 increments of thick- 
ness Ar®. Two grid stations were located at the interface, one corresponding 
to the 3-phase and the other to the a-phase. The rate of travel of the nth 
grid point is related to the velocity of the interface by 

^ ^ d(^/2) (7 

dt ^/2 dt 


for the nth point in the 3-phase, or by 
drn ^ /L/2 - rn \d(^/2) 

dt Il/ 2 - C/2/ dt 


for the nth point in the a-phase. The rate of change of concentration at any 
internal station can be expressed as 

*^*^n = ^*^n ^^n + ^^n (c 

dt 3rjj dt 9t 

where SC^/Bt is given by equation (1). By combining equations (7) and (9), 
an expression for the variation of concentration with time at the nth internal 
station in the 3-phase can be written as follows: 

^*^n _ ^n ^*^n d(^/2) + ^*^n C1C 

dt (^/2) 3rn dt 3t 


In an analogous fashion. 


*^^n _ /l/ 2 - rji \9Cn d(^/2) + ^^n 
dt ~ I L/2 - 5/2/3rn dt 3t 


gives the change in concentration with time at the nth internal station in the 
a-phase . 


These equations can be normalized by the following change of variables: 


Rn = 


T _ ^WxL 
(L/2)2 
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By using these variables, equations (10), (11), (1), and (6) were rewritten as 
3 -phase: 

d(^ ^ ^ ^ d(^/L) ^ ^ 


dT K/L 9Rn dT 3T 


( 12 ) 


ct-phase: 


where 


d*^n - / l - Rg \ ^^n d(C/L) + ^*^n 
dX \1 - ^/L/3Rn dX 3X 


3X 3Rn 3Rn 3Rn 


(13) 


d(C/L) . _l_ 
dT 63a 




(14) 


(15) 


Xhe symbol 63(;;j denotes the difference between the solubility limits of 6- and 
ct-phases . 


Equations (12) to (14) were rewritten in finitfe-difference notation (expli- 
cit form, second-order central difference) as follows: 


3-phase: 

Cn"^ - Cn _ ^ 

AX C/L dX y 2 Ar3 / 9T 

a-phase: 

- Cn - h - Rn ^Cn-Hl - C^-l\d(C/L) ^ ^ 
AX \^1 - K/lI\ 2 AR« / dX 3X 


where 

9Cn . nJlm (cL^ ~ cil\ , (cL^ - 2C^ + C^.A 

3X 2 AR j \ (ar)2 j 


+ 




nJ 


^n +1 ~ ^n-1 
2 AR 


^ Cn-hl - Cn-1 
( 2 AR y 


( 16 ) 


(17) 


( 18 ) 


Xhe finite-difference expression (second-order forward and reverse difference) 
for the interface-flux-balance equation (15) is 
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(19) 


AT 


1 

<S6a 


- %a 


/-Ci+3 + 4 Ci+2 

T > 

- 31^6 

V 2 ARpt 

! 


f^I-2 ■ ^^1-1 * 


2 AR® 


J 


The finite-difference expressions for the boundary conditions are given by 

j+1 J i /J J' 

= 2(m + 1)5i(?2_l^\ (R = 0) (20) 

\(AR3)^ 


pJ' 

- Cn 

AT 


= 25n 


- Cn 


(AR“)' 


(R = 1) 


( 21 ) 


To perform the finite-difference calculations, the grid spacing AR must 
be chosen small enough to give convergent solutions of acceptable accuracy. 
Also, because explicit formulation was used in the development of the equa- 
tions, the time increment AT must be selected according to the stability 
requirement (similar to eq. (5)), 

AT^x ^ 0.25(ARmin)2 ^^2) 

where AT^ay is the maximum time increment and ARjnin is minimum grid 
spacing. 

A computer code was developed to solve equations (16) to (22) for any two 
phase binary alloy system. The FORTRAN IV code, documentation describing the 
code, and a sample case illustrating typical input and output information are 
available from COSMIC (ref. 9). The special features of the analysis included 
in this code are discussed in a subsequent section. 


Three-Phase Systems 


Diffusion between pure metals which form a two-phase binary alloy system 
with an intermediate phase will result in a composition variation between the 
two metals of the form shown in figure 1(c). The concentrations at the 3 y 
interface correspond to the equilibrium solubility limits of the 3~phase (Cgy) 
and the y-Pbase (Cyg) at the temperature of interest. The concentrations at 
the ycx interface correspond to the equilibrium solubility limits of the y-phase 
(CyQj) and the a-phase (C^y). 


Diffusion in this system can be described by Pick's second law for each 
phase and two interface-flux-balance equations. The interface-flux-balance 
equations are of the form 




d(^l/2) 

dt 



r<*/2 


■ 


r=^?/2 


(23) 
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at the yot interface, where ^i/2 is the location of the By interface with 
superscripts + and - denoting the y and B sides of the interface, respec- 
tively, and ^ 2/2 is the location of the yot interface with superscripts + and 
- denoting the a and y sides of the interface, respectively. The initial 
and boundary conditions to be satisfied are 


C = C (0 1 r < !l/2; t = 0) 


C = Co (i/2 < r i L/2; t = 0) 


= 0 (r = 0 and L/2; t ^ 0) 

9r 

All interface concentrations are assumed equal to the equilibrium solubility 
limits of the phases : 

By interface: 

C = C3y = ^?/2) 

C = Cyg (r = K*/2) 

ya interface: 

C - Cyot (r = ^2^2) 


C = Ccty (r = K 2 / 2 ) 

In general, the thicknesses of the B-, y-, and ot-phases change with diffu- 
sion time. To account for these dimensional changes, the variable-grid finite- 
difference technique previously discussed for the two-phase analysis was used. 
The number of grid points in each phase was held constant and the grid sizes 
ArP, Ar’Y, and Ar^ changed with time. The rate of change of concentration 
at any internal station is given by equation (9) where the velocity of the nth 
grid point dr^j/dt in each phase is given by 

^*^n - *^n /2) ^25 

dt C 1/2 dt 


for the nth point in the B-phase, by 

dr„ ^ d(Ci/2) ^ |"d(C 2 / 2 ) _ d{K^/2)' 

dt ~ dt "L dt 
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for the nth point in the yP^^sse where 


K 


n 


^n - ^1/2 
^2/2 - Cl/2 


and by 

dr„ / L/2 - r„ \d(^2/2) 

dt " (l/ 2 - ^2/2/ dt 


( 27 ) 


for the nth point in the a-phase. The equations for the rate of change of con- 
centration at any internal station in the three phases are obtained by substi- 
tuting equations (25), (26), and (27) into equation (9): 


3 -phase : 


dC 


n _ 


dt 


*^n d(^-|/2) 9C|^ ^ ^('n 

CT/2 ~dt ^ aT" 


( 28 ) 


y-phase : 


, /d(Cl/2) ^ . 


dt 


dt 


d(C2/2) _ d(^i/2)‘ 


dt 


dt 


3rn at 


(29) 


a-phase : 


dCn /L/2 - r„ \d(^2/2) ac^ ^ ^ 

dt \L/2 - 52/2/ dt ar„ at 


(30) 


Normalizing the distance, time, and diffusion coefficient in the same manner as 
was done previously for the two-phase analysis reduces equations (28) to (30) to 


3 -phase : 

dCn _ Rn d(5i/L) ac^ , aCn 
dT 5i/L dT 3Rn 3T 

y-phase : 

dCn _ jd(5i/L) ^ ^ rd(52/L) . d(5i/L)l'l aCn ^ aq 
dT " ^ dT "L dT * dT J J 3Rn 9T 

a-phase : 

d^ - / 1 - Rn V(C2/L) ^ ^ 

dT - 52/L/ dT 3Rn ai 


(31) 


(32) 


(33) 


where aC(i/3T is given by equation (l^l). Likewise, the interface-flux-balance 
equations (23) and (24) can be rewritten as 
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By interface: 
dT 


T 


ya interface: 




dR L_^ + 


R<T/L 


R=5l/L 


E^2^L 


where 


<SBy = ^By “ S'B 


Sy(x = Cyot - Cq[y 

Equations (31) to (35) can be expressed in finite-difference notation in 
a manner analogous to that discussed for two-phase systems: 


3 -phase : 


^ ~ ^n - ^n d(^i /L) /cp+l ~ ^n-l\ + ^^n 


Y-phase : 


5l/L dT 


2 AR^ 


- ci _ fd(^i/L) ^ [d(^2/L) _ d(^i/L)]'\ Zc^^i - Ci_\ ^ 30^ 
AT dT " L dT dT _J 2 ARY / 3^ 


a-phase : 


pj+’l pj 


1 - Rnya2/L)(C^+1 - C^-l\ , 3C„ 

1 - C2/W dT I 2 AR« / 9T 


where 3Cn/3T is given by equation (18). The finite-difference expressions for 
the interface-flux-balance equations (3^) and (35) are 

By interface: 

(Ci/L)'^'^^ - (Ci/L) -^ ^ 1 g g(^ -cj-|+3 + ^^i^+2 - 3C^3 ^ 

At ^ \ 2 ArY ^ 

- 5e-/ ^Il-2 - ^^Ii-1 ^ (39) 

2 AR3 } 
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ya interface: 


(^2/L)'^'^^ - _ 1 


AT 


^yaL 


Dayf"^^2+3 + '*Ci2+2 - 3Ca^) 


2 AR“ 


- Dya 


^Cl2-2 - ^Ci2-i 
^ 2 Ar'Y 


3CyO^ 


( 40 ) 


The finite-difference expressions for the boundaries are given by equations (20) 
and (21 ) , 


The grid spacing AR was chosen small enough to give convergent solutions 
of acceptable accuracy, and the time increment was chosen according to the same 
stability requirement that was used for the two-phase system (eq. (22)). 

A computer code was developed to solve equations (18), (20) to (22), and 
(36) to (40). The FORTRAN IV code, documentation describing the code, and a 
sample case illustrating typical input and output information are available 
from COSMIC (ref. 10). The special features of the analysis included in this 
code are discussed in the following section. 


SPECIAL FEATURES OF ANALYSES 
Single-Phase Analysis 

Volume change .- To account for volume changes which occur during diffusion 
due to differences in the molal volumes of the pure metals, the concentration 
was expressed in terms of the volume fraction of one of the diffusing components. 
All concentrations were converted into volume fraction for the calculations and 
then converted back to atomic fraction once the desired solution had been 
obtained . 

Logarithmic interpolation of D .- Logarithmic interpolation of the 
concentration-dependent diffusion coefficient was employed in the analysis 
because of the general trend for log D to vary linearly with C for several 
binary alloy systems (ref. 13). 

Initial composition profile .- Two optional initial conditions were treated: 
(1) no prior diffusion between the A-rich and B-rich regions, and (2) prior 
diffusion resulting in a concentration gradient across the interface region. 

Most analyses treat the case of no prior diffusion where both regions of the 
couple are uniform in composition with a sharp change in composition at the 
interface. The second option, however, allows an experimentally determined 
concentration profile to be treated as the initial condition to calculate how 
much more diffusion will occur with additional exposure at elevated temperature. 
This feature was added because in the consolidation of many metal-matrix compos- 
ites, some diffusion occurs between fibers (or particles) and the matrix. 

Finite-difference grid size .- For long diffusion times, one or more grid 
changes were made to reduce the total number of grid stations to a minimum (20) 


12 


in an effort to minimize computer time. The concentrations corresponding to 
each new grid point were determined by interpolation between points on the old 
grid. 


Two-Phase Analysis 

Differences in molal volumes, concentration dependence of the diffusion 
coefficient, and the presence of an initial concentration profile in the sample 
were accounted for by using the procedures previously discussed for single- 
phase systems. 

Finite-difference grid size .- To improve solution accuracy, grid stations 
were concentrated near the interface during the initial stages of diffusion 
when the diffusion zone was small relative to the specimen dimensions. When 
the concentrations at the outer bounds of the zone being analyzed began to 
change, the zone size was enlarged. In general, two or three zone expansions 
were performed before the bounds of the a- and B-phases were reached. In all 
cases, the number of grid stations in each phase was held constant. After each 
expansion of the zone, the concentrations at the new grid locations were deter- 
mined by interpolation from the concentration profile (present in the specimen) 
just prior to expanding the zone. 

Interface-location criteria .- For the two-phase system, solution accuracy 
was improved by (1) employing an iterative mass-conservation procedure to locate 
the aB interface rather than relying on the interface-flux-balance relation- 
ship and (2) using logarithmic average diffusion coefficients at the interface 
rather than the diffusion coefficients corresponding to the concentrations at 
the interface. Both of these changes compensate for numerical approximation 
errors inherent in the procedure used to calculate the interface concentration 
gradients . 

Logarithmic average interface diffusion coefficients were calculated by 
the expressions 

^aB = hoL = l/^I-I^I 

where and are the normalized diffusion coefficients used in 

equation (19). 

The iterative mass-conservation procedure consisted of solving the interface- 
flux-balance equation (eq. (19)) to obtain an estimate of the location of the 
aB interface and adjusting this location in the positive or negative direction 
to force the system to conserve. This refined interface location was used to 
calculate new grid spacings AR*^ and AR^ which in turn were used to calculate 
new grid-point locations. The concentrations at the new grid points were assumed 
equal to the concentrations at the old grid points before shifting the interface. 
The mass-conservation procedure was repeated to obtain a further refinement in 
the interface location. In general, only two or three iterations were required 
to establish convergence . 
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Phase -termination criteria.- For all B-phase concentrations less than 


^Ba ^^min stability criterion (eq. (22)) was taken as Afl®, 

and the B-phase concentrations were not computed from this time onward. In an 
analogous fashion, the a-phase concentrations were not computed if all the 
a-phase concentrations were greater than - 10"^. If both conditions were 

met, or if the B-phase thickness was 0.01 of its initial thickness, no further 
calculations were performed. 


Three-Phase Analysis 


To account for differences in molal volumes, concentration dependence of 
the diffusion coefficient, and the presence of an initial concentration profile 
in the three-phase system, the procedures previously discussed for single-phase 
systems were employed. As previously discussed for the two-phase analysis, grid 
stations were concentrated near the interface during early stages of diffusion, 
and the zone of analysis was expanded with diffusion time. 


Initiation and growth of intermediate phase .- Initiation of the intermedi- 
ate phase (y) was treated by assuming that y was initially present in the 
sample in the form of a thin (0.05 ym) layer between the a- and B-phases. The 
concentration gradient in y was assumed to vary linearly between the solubil- 
ity limits Cyg and Cycc If was less than D® or D^, the y-phase 

remained thin 'for most or all of the homogenization process depending on the 
value of the average composition C and the solubilities of the a- and B-phases. 
Under these conditions, a large amount of computation time was required to per- 
form the calculations because the grid size in the y-phase Ar"'^ remained small 
and dictated a very small At (eq. (22)). To overcome this problem, the 
inter face-flux -balance equations were solved to determine the thickness change 
of each phase with time, but compositional changes in the y-phase were not com- 
puted and a linear gradient was assumed. Because y-phase concentrations were 
not calculated, the time increment for the finite-difference analysis was deter- 
mined by the relatively larger grid spacings in the a- or B-phases. This proce- 
dure resulted in significant savings in computation time. If the y-phase thick- 
ness grew larger than a small preset value, concentrations in the y-phase were 
calculated . 


To prevent oscillations in the thickness of the y-phase, the time increment 
for the diffusion calculations was selected small .enough to limit the thickness 
change of the y-phase to no more than 2 percent of its previous value during 
any one time irtcrement. This was accomplished by calculating At for the a- 
and B-phases and using the smaller of these to calculate the change in thickness 
of the y-phase. If the change was more than 2 percent. At was reduced to give 
a 2-percent change in the thickness of the y-phase. 

Phase-termination criteria .- If any phase thickness was less than 1 percent 
of the total zone thickness, or if all the concentrations within a given phase 
were within 10~”, the concentrations in that phase were assumed to remain con- 
stant during subsequent time increments. Also, the grid spacing in this phase 
was not used to compute the time increment . If the thickness of any phase was 


less than 0.001 percent of the total diffusion-zone thickness, this phase thick- 
ness vfas assumed to be equal to zero. 


RESULTS AND DISCUSSION 
Mass Conservation 

As previously discussed, an iterative mass-conservation procedure was 
used in the two-phase analysis instead of the conventional interface-flux- 
balance procedure. Calculations were performed to assess the improvement in 
solution accuracy resulting from this change. Typical results are shown in 
figure 2. The conditions selected for this assessment were: cylindrical 
geometry, S./2 = 20 ym, L/2 = 140 ym, C* = 1.0, Cga = 0.85, CQ[g = 0.15, 

Cq = 0, and D°^ = D® = 10“^^ m^/s. Solutions were calculated for initial cx-phase 
grid sizes of 0.5, 0.86, 3, and 6 ym and an initial 8-phase grid size of 1 ym for 
all cases. Figure 2(b) compares normalized interface locations at t = 90 000 s 
calculated by the two interface-location criteria as a function of the initial 
grid size in the ot-phase. Use of mass conservation resulted in a faster rate of 
convergence as the a-phase grid size was decreased. Very small initial a-phase 
grid sizes were required to predict the same interface location by both methods. 
If the flux -balance criterion were used, very small grid sizes are required to 
give conservation of mass; see figure 2(a). These results suggest that the mass- 
conservation criterion improves solution accuracy and significantly reduces com- 
putation time because a larger grid size can be used. 

A mass-conservation criterion cannot readily be used in the three-phase 
analysis because of the additional complexity associated with the presence of 
two interfaces. If the system deviates from conservation, no simple way was 
found to determine how much each interface should be shifted to force the system 
to conserve. 


Intermediate-Phase Growth 

The procedure used to treat the initiation and growth of an intermediate 
phase in the three-phase analysis was previously discussed. If the intermediate 
phase was thin, the three-phase problem was essentially reduced to a two-phase 
problem by not calculating the concentrations at the finite-difference grid sta- 
tions located in the intermediate phase. However, the time increment At for 
the finite-difference calculation was selected small enough that the growth of 
the intermediate phase was stable. The At was chosen so that the thickness 
of the intermediate phase did not change by more than 2 percent during any time 
increment. The effectiveness of this procedure is illustrated by the results 
shown in figure 3. Compared in this figure are Y-phase thicknesses determined 
with (1) the three-phase calculation, (2) the two-phase calculation using a 
±2-percent-change criterion to select At, and (3) a two-phase calculation 
using a change criterion of -50 to +100 percent to select At. The three-phase 
calculations and the two-phase calculations with a ±2-percent-change criterion 
converge after approximately 60 s of diffusion time. The two-phase calculation 
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with the change criterion of -50 to +100 percent was unstable and caused the 
Y-phase thickness to oscillate about the three-phase results. 

For the diffusion conditions selected for this illustration, the computer 
time required to calculate the solution for 60 s of diffusion time was more than 
a factor of 100 larger for the three-phase calculation than for the two-phase 
calculation using the ±2-percent-change criterion. Although this time differ- 
ence was somewhat dependent on the solubilities and diffusion coefficients of 
each phase, the difference was always large for those cases considered where 
d'Y < D« or d3. 


Example Cases and Comparisons With Other Solutions 

Single-phase analysis .- A comparison of results calculated by the present 
finite-difference solution and an iterative solution technique (ref. 4) also 
capable of treating concentration-dependent D is shown in figure 4. In this 
calculation, D varied linearly on a logarithmic scale with concentration. 

The D for pure A was assumed to be an order of magnitude larger than the D 
for pure B. The exposure time was 5 hr and the geometry of the specimen was 
planar. The iterative solution technique of Unnam and Houska (ref. 4) is appli- 
cable only for planar geometry. The excellent agreement between the two solu- 
tion techniques is typical of that found for all other exposure conditions con- 
sidered. These results provide a check for the finite-difference analysis 
developed in this investigation. 

Two-phase analysis .- Results of this analysis were compared with those of 
a closed-form iterative technique reported by Unnam and Houska (ref. 4). Com- 
parable solutions were obtained for a wide range of input parameters for planar 
interfaces (solution of ref. 4 is valid only for planar interfaces). Several 
different phase thicknesses and D variations were considered in order to yield 
information regarding solution stability and accuracy. No significant differ- 
ences were found between the results obtained from the two solutions. The inter- 
face positions calculated by the finite-difference program were always within 
5 percent of the positions calculated by the iterative solution. 

A comparison was also made with experimentally determined diffusion data 
on W-Ni laminates reported by Tanzilli and Heckel (ref. 2). They prepared mul- 
tilayer W-Ni couples with two mean compositions, C = 0.152 and C = 0.121 
atomic fraction of W. These couples were diffused at 1480 and 1429 K, and then 
they were sectioned perpendicular to the layers and metallographically polished. 
The extent of diffusion was determined by measuring the average thickness of the 
W-rich layers after each exposure. Their results are presented in figure 5 
where the normalized thickness of the B-layer, ^/l, is plotted as a function 
of a dimensionless time factor, where ^ is 6-phase (W-rich) 

thickness at time t, I is the initial thickness of this layer, and IP is 
the concentration-independent diffusion coefficient in the ot-phase (Ni-rich). 

Also plotted on this graph are results obtained from (1) a constant-D finite- 
difference solution developed by Tanzilli and Heckel (ref. 2), (2) the present 
variable-D finite-difference solution with the planar geometry option, and 
(3) the iterative solution developed by Unnam and Houska (ref. 4). For cases (2) 
and (3), the diffusion-coefficient data of Walsh and Donachie (ref. 14) were 
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used for the a-phase region and the data of Vasilev and Chernomorchenko (ref. 15) 
were used for the 3 -phase region. The same initial and boundary conditions were 
used in all these solutions. The analytical results are shown only for the sam- 
ple with C = 0.152, as the curves were the same for C = 0.152 or 0.121 until 
C/fi- dropped below about 0.3. Also, no experimental data points were available 
below ^/i = 0.3 for the C = 0.121 sample. 

The interface positions calculated by all three solutions agree. The curve 
obtained from the present analysis best fits the experimental data points. 

Three-phase analysis .- To study the effect of large variations on the 

diffusion process, calculations were performed for EP^/D^ ratios of 1, 100, and 
0.01 holding and constant at 10“^^ m^/s. Average concentration of the 

couple was 0.40 atomic fraction of B and the solubility limits were C = 1.0, 

Cgy = 0.99, Cyg = 0.51, Cycj = 0.49, C<xY " 0.20, and Cq = 0. The small solu- 

bility range of the B-phase, 0.01, and Y-phase, 0.02, mean that any differences 
in the diffusion kinetics are predominantly due to differences in diffusion 
in the a-phase. The normalized thickness of the 3- and Y-phases are plotted 

as a function of \Jt in figure 6(a). The curves for eP = Dp = D ' provide a 

basis for comparison of the other two curves. If EP^ is small (D<^/d 3 = 0.01), 
the rate of decrease of the 3 -phase is slower and the rate of growth of the 
Y-phase is faster than the corresponding curves for z 1 . If iP' is 

high (EP^/D^ z 100), a rapid loss of 3-phase thickness and slow initial growth 
of the Y“Phase occur. At approximately Vt = 50, the rate of decrease of the 
3-phase is greatly reduced and the rate of growth of the Y“Phsse increases. 

The corresponding concentration profiles (fig. 6(b)) show that these changes 
are due to the small concentration gradient in the a-phase for the high D® 
case. A small gradient in the a-phase implies a small flux into the Y-phase 
which increases the growth rate of the Y-Phase and reduces the rate of loss 
of the 3 -phase. 


CONCLUDING REMARKS 

Finite-difference solutions of the one -dimensional transient diffusion 
equations describing diffusion in single-, two-, and three-phase binary alloy 
systems were developed. These solutions are applicable for any diffusion-zone 
size and for planar, cylindrical, or spherical geometries. Each analysis also 
treats concentration dependence of the diffusion coefficient. General computer 
codes were developed to perform each analysis. These computer codes, documenta- 
tion describing the codes, and sample cases illustrating typical input and out- 
put are available from COSMIC. Special features were included in the computer 
codes to account for differences in molal volumes, initiation and growth of an 
intermediate phase (three-phase system), disappearance of a phase (two- and 
three-phase systems), and the presence of an initial composition profile in the 
specimen. Each computer code was optimized to achieve good accuracy while mini- 
mizing computation time. The use of a mass-conservation criterion to locate the 
interface in the two-phase system rather than the conventional flux-balance 
criterion resulted in improvement in solution accuracy and reduction in computa- 
tion time. In the three-phase analysis, computation time was reduced by not per- 
forming the concentration calculations in the intermediate phase when that phase 
was thin. In all three computer codes, the grid spacings were increased during 
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the later stages of diffusion when the concentration gradients were small. The 
increased grid spacing allowed increased time increments and thus minimized com- 
putation times. In the single-phase analysis, this was accomplished by reducing 
the number of grid stations with diffusion time. In the two- and three-phase 
analyses , a given number of grid stations were initially concentrated in a nar- 
row zone near the interface. For longer diffusion times, this zone was expanded 
until the specimen boundaries were reached. 

The results of the single- and two-phase analyses were in good agreement 
with solutions or data reported in the literature for those cases where a direct 
comparison could be made. A direct comparison of the three-phase results was 
not possible because neither experimental nor accurate analytical data were 
available from the literature. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
January 27, 1978 
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(a) Single phase - complete solid solubility. 



(b) Two phase - limited solid solubility. 



(c) Three phase - intermediate-phase formation. 

Figure 1.- Phase diagrams and example concentration profiles produced by dif 
fusion between initially pure A and pure B at temperature 6j) in single- 
two-, and three-phase binary alloy systems. 


Interface-location 

criterion 


Interface flux balance 

Mass conservation 



Initial grid size In a-phase, |im 


(a) Mass conservation. 



Initial grid size In a-phase, pirn 


(b) Convergence of interface location. 

Figure 2.- Effect of interface-location criteria on (a) mass conservation 
and (b) convergence of interface location for different initial grid 
sizes in a-phase. Cylindrical geometry. 
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Figure 4.- Comparison of concentration profiles in single-phase system calcu- 
lated by the present analysis and an iterative solution reported by Unnam 
and Houska (ref. 4). Planar geometry; = 10D®; t = 5 hr. 
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Figure 5.- Comparison of calculated and experimental interface migration data 
for planar W-Ni couples. C is given in atomic fraction of W. 



(a) Normalized thickness of 3 - and y-phases. 



0 20 40 60 80 

Distance, mhi 

(b) Concentration profile at \ft = 50. 


Figure 6.- Effect of large iP' variations on (a) the change in normalized 
thickness of 6- and Y-phases and (b) the concentration profiles for a 
cylindrical geometry specimen with C = 0.40 atomic fraction of B. 
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